EPJ manuscript No. 

(will be inserted by the editor) 



3 ' Plasticity and dynamical heterogeneity in driven glassy 
<^ ■ materials 

c : 

■ Michel Tsamados 
■ 

^\ ' Universite de Lyon; Univ. Lyon I, Laboratoire de Physique de la Matiere Condensee et Nanostructures; CNRS, 
CSJ ! UMR 5586, 69622 Villeurbanne, France 

, Received: date / Revised version: date 

. 

, Abstract. Many amorphous glassy materials exhibit complex spatio-temporal mechanical 

■ response and rheology, characterized by an intermittent stress-strain response and a fluctu- 
' ating velocity profile. Under quasistatic and athermal deformation protocols this heteroge- 

^ , neous plastic flow was shown to be composed of plastic events of various sizes. In this paper, 

1 . through numerical study of a 2D Lennard- Jones amorphous solid, we generalize the study of 
' O ' the heterogeneous dynamics of glassy materials to the finite shear-rate (7 0) and temper- 
Ch I ature case (T 7^ 0). The global mechanical response obtained through the use of Molecular 
O . Dynamics is shown to converge in the limit 7 — >■ to the quasistatic limit obtained with 

I ^ I ' an energy minimization protocol. The detailed analysis of the plastic deformation at differ- 
ent shear rates shows that the glass follows different fiow regimes. At sufficiently low shear 

^-H ■ rates the mechanical response reaches a shear-rate independent regime that exhibits all the 

^ ' characteristics of the quasistatic response (finite size effects, yield stress...). At intermediate 

, shear rates the rheological properties are determined by the externally applied shear-rate. 

■ Finally at higher shear the system reaches a shear-rate independent homogeneous regime. 
\l ' The existence of these three regimes is also confirmed by the detailed analysis of the atomic 

motion. The computation of the four-point correlation function shows that the transition 

. from the shear-rate dominated to the quasistatic regime is accompanied by the growth of a 

' dynamical cooperativity length scale ^ that is shown to diverge with shear rate as ^ oc 7"", 

, with 1/ ~ 0.2 — 0.3. This scaling is compared with the prediction of a simple model that 

7—i ■ assumes the diffusive propagation of plastic events. 

> ■ 

. ^ , PACS. PACS-key discribing text of that key - PACS-key discribing text of that key 

H : 

. . , 1 Introduction features can be seen as the macroscopic manifesta- 
tions of the strongly heterogeneous underlying de- 
When subjected to slow driving many systems ex- formation processes that take place in the driven 
hibit an mtermittent response with the appearance gj^^g^gg ^^^^ je^^ ^^.^^-^ velocity lo- 
of discrete and impulsive events spanning a broad calization observed experimentally in alloys, metal- 
range of sizes. Such a scale-mvariant behavior is ^-^ gj^^^ggg ^ polymers ^ > granular media [U], 
generally observed m driven nonlinear, dynamical ^^^^g ^^^^^-^^^ |4l[T2l[T3l[T4| as well as in nu- 
systems and examples of such crackling signals are ^^^^^^ simulations, both of model systems such as 
ubiqmtous m nature (for a review see [T]). Lennard- Jones glasses [T51[TC1|T71[T5] , as in more re- 

Among such systems varioi^ glassy materials ^listic simulations 
(granular media [2j|3], foams [Hlo], emulsions [HIS], 

micelles [7], metallic glasses [!]■•■) were shown to Experimental studies on disordered materials (foams, 
exhibit, at small scales, small shear rates and for granulars), far below the glass transition temper- 
sufficiently low temperatures, such an intermittent ature P51[^[T^[^[^[?7[|28[ H1[T ^ . have associated 
signature in their stress-strain mechanical response the dissipative events observed in the stress strain 
along with complex spatio-temporally fluctuating mechanical response (the stress drops) to the exis- 
velocity profiles. Such intermittent spatio-temporal fence of a collective behavior of localized rearrange- 
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ments leading to a strongly heterogeneous mechan- 
ical response as shown in the existence of shear 
bands in the macroscopic plasticity of such sys- 
tems. Along these experimental evidences, in the 
last five to ten years, extensive numerical simula- 
tions of quasistatically sheared model glassy sys- 
tems have confirmed this picture of a cascade mech- 
anism but have also generated some debates first as 
to the validity of the potential energy minimization 
(PEM) methods to represent the physical reality 
of slowly driven systems and second as to the lo- 
calized nature j29j of the elementary rearranging 
plastic events associated to the macroscopic stress 
and energy releases of figure [H 

On a theoretical level various current 'mean field' 
models of the rheology of glasses predict reasonably 
well the macroscopic mechanical properties of these 
materials [5U1I5T] . However to reproduce the spa- 
tial heterogeneity of their response a novel class of 
spatially resolved models inspired from geophysics 
[32] have depicted the plastic flow as a result of the 
collective organization of interacting plastic rear- 
rangements in an otherwise elastic medium. This 
type of model allows therefore a mesoscopic ap- 
proach of the appearance of spatio-temporal het- 
erogeneities based on elementary bricks, the plas- 
tic rearrangements. Nevertheless to produce a mul- 
tiscale description of the mechanical response of 
glassy materials necessitates a microscopic (i) iden- 
tification of these elementary plastically rearrang- 
ing zones |18 p 33 ( [34 ] . (ii) a comprehension of their 
local elasto-plastic properties [3^1 and yield crite- 
ria [3 6) and (iii) a study of their interaction and 
propagation. 

Basing our study on the detailed numerical anal- 



acterize the spatio-temporal dynamics of the driven 
glasses we first compute various two-time observ- 
ables (mean-square displacement, van Hove func- 
tion, intermediate scattering function...) which in- 
form on the typical relaxation times and their link 
with the external shear rate, we then proceed to a 
general description of the dynamical heterogeneity 
of the driven systems and prove the existence of a 
diverging dynamical cooperative length scale as the 
shear rate 7 tends to zero. We also provide pictorial 
evidence of these dynamical structures for different 
shear rates. Finally we discuss our results in light 
of recent advances and compare our observations to 
a phenomenological mesoscopic yield stress model. 



2 Rheological response of the system 

2.1 Rheological and mechanical characteristics. 

We do not reproduce here the details of the qua- 
sistatic simulation procedures that we use in this 
paper as they were presented already in details in 
[37ll38l[39lfT8] . Here we extend the results obtained 
on the quasistatically and athermally sheared two 
dimensional polydisperse glasses presented in our 
earlier work to finite shear rates by the use of molec- 
ular dynamics while the temperature is maintained 
well below the glass transition temperature. As in 
|14| we impose the transverse temperature {Ty = 5 • 
10~^) through the use of a simple velocity rescaling 
thermostat on the transverse component of the par- 
ticle velocities. For the applied shear an equivalent 
amount of runs were produced under two bound- 
ary condition protocols, namely rigid walls bound- 
ary conditions (simplified notation RWBCs) at y 



ysis of a sheared two dimensional polydisperse Lennardi^/2, where H is the height of the sample, and 



Jones glass [571I551[551|181I36[I55] we discuss mainly 
in this paper point (iii) with an emphasis on the 
influence of the shear rate on the cooperative na- 
ture of the plastic response. At a microscopic level a 
plastic rearrangement redistributes its shear stress 
through a quadrupolar long-range elastic propaga- 
tor [40ll4T| hence favoring subsequent plastic rear- 
rangements in its vicinity and within some pre- 
scribed directions. An estimate of the dynamical 
coopcrativity length scale resulting from this cas- 
cading process can be computed with tools such as 
the four-point correlation function borrowed from 
the literature on dynamical heterogeneities in su- 
percooled liquids near their glass transition [42] or 
in granular media near their jamming transition 

In the present paper we provide an extensive 



Lees-Edwards boundary conditions (noted LEBCs). 
In total we have analyzed here 24 samples of (625 
particles, Lx=25.9938, Ly=25.9938), 8 (2500,51.9875,51.9875) 
and 8 (10000,103.975,103.975) aU corresponding to 
a density p = 0.925. We have sheared at finite shear 
rates also one larger sample (40000,206.950,206.950) 
that was not sheared with the quasi-static protocol. 
For each simulation, we collect data over four strain 
units (e = 4) and store all the positions of the par- 
ticles at a regular strain interval of 6e = 10~^. In 
figure [T] we have reported the mechanical response 
of the smallest sample for RWBCs. We observe the 
characteristic fiow behavior associated with glassy 
materials with a convergence of the response to the 
quasistatic limit as the shear rate is progressively 
reduced and a global non linear flow curve of the 
Herschel-Bulkley type : t = Ty + ci'y'^, where ry is 



study of the dynamics and plasticity in a two-dimensioiikl; yield stress, and where one can define the vis- 
model Lennard- Jones glass. To quantitatively char- cous stress Ty = t — Ty. In figure [2] we draw these 
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flow curves for both boundary conditions and for all 
system sizes. The stress-strain rate curves of figure 
[2] are obtained by averaging the macroscopic stress 
values of the sheared glass obtained in figures [T] for 
strains larger than e > 25%, i.e. deep in the plastic 
flowing regime once the stationary flowing regime 
is established. Indeed we have checked that a linear 
velocity profile is established in the different sam- 
ples for typical strains of the order of e ^ 2.5%. 

2.2 Convergence to the quasistatic limit. 

Importantly the results presented above bridge the 
gap between the two types of approaches used in the 
literature, namely quasistatic energy minimization 
protocols and finite shear rates molecular dynamics 
methods, and resolve the controversy relative to the 
validity of the quasistatic protocols. Indeed as is ev- 
ident from the mechanical response shown in figures 
[T] the quasistatic stress-strain curve appears as the 
limiting curve of the finite shear rates procedures. 
The superposition of the quasistatic response with 
the 7 = 10~^ shear rate response is in fact almost 
perfect in the early parts of the curves before small 
differences are amplified irremediably. One sees for 
example on the right hand side of figure [1] how even 
very small relaxations at around e ~ 3.2 — 3.4% in 
the quasistatic response (black thick line) are also 
visible in the lowest shear rate curves (yellow and 
magenta). Note that these small features of the me- 
chanical response of the glass would not be visible 
if the temperature was higher and therefore induc- 
ing a noisier stress signal. The good convergence to 
the quasistatic protocol is also apparent in figures [5] 
were one can indeed observe that the values of the 
lowest finite shear rates are in good agreement with 
the 7 = shear rate method (under both boundary 
conditions <Jxy{j < 10~^) ^ 0.4). again confirming 
the physical relevance of the quasistatic method. 

2.3 Decomposing the plastic events in 
elementary units. 

In an important scries of papers |33[I34I it was re- 
cently shown by Lemaitre and Maloney that the re- 
laxation in mechanically driven glasses occurs through 
the formation of cascades of quadrupolar elemen- 
tary units. In the quasistatic protocol to decompose 
the cascade in its subunits one needs to study in 
detail the evolution of the positions of the particles 
during the minimization procedure, with the inher- 
ent limitation of the minimization algorithms that 
one can not associate a time scale to the succes- 
sive elementary rearranging units. This limitation 



is automatically overcome in the molecular dynam- 
ics simulation where time appears explicitly in the 
algorithm and where one can follow the evolution of 
the cascade in time. As is apparent at finite shear 
rates in figures [1] the plastic relaxation is not in- 
stantaneous during a stress drop and this typical 
lifetime of the plastic events is transposed in the 
stress-strain mechanical response in the downward 
portions with the slopes that increase with increas- 
ing strain rate. To understand intuitively the mech- 
anisms involved in the mechanical response of the 
sheared glasses and the different time scales that are 
relevant it is highly instructive to visualize movies 
of the instantaneous non-affine displacement field 
during the deformation. The two panels of figures 
[4] and [5] show snapshots of such movies taken at 
regularly spaced strain intervals corresponding to 
the red symbols on figure [3l As can be seen in this 
last figure, for the slowest strain rate 7 = 10~^ the 
strain interval between snapshots is Se = 2.5 x 10"'^ 
(red triangles) and for the fastest strain rate 7 = 
10^^ it is Se = 1.25 x 10^^ (red circles). When look- 
ing at the two panels one must bare in mind that 
while in the slowly sheared case the total strain ap- 
plied between the first figure and the last is less than 
1% in the fast case it is more than 5%. Also the non- 
affine displacement field represented on each snap- 
shots correspond to the displacement during a time 
interval of St = ILJU and therefore the associated 
strain S"f ~ x St is one hundred times larger for 
the fast shear than for the slow shear. If one as- 
sumes that the density of weak triggering zones of 
plasticity is homogeneously distributed through the 
sample with a shear rate independent density per 
unit strain f2 then one expects also to see 100 times 
more such nuclei of plasticity in the fast case. This 
explains the visual impression that there exists a 
higher density of local plastic displacements (black 
arrows) on panel [5] than on panel |4] where only a 
few local plastic rearrangements are observed. An- 
other difference between the slowly and fast driven 
regimes is that while in the fast case the nuclei 
do not seem to merge or percolate or evolve in a 
correlated manner in the slowly driven case the lo- 
cal quadrupoles show a cooperative and correlated 
avalanche dynamics. This is particularly visible on 
the fourth panel of figure|4]where one can see about 
^ 5 such elementary plastic units forming a L-shape 
with four events horizontally aligned. Interestingly 
the typical distance between quadrupoles on this 
figure is about ^ ~ 20LJU which is reminiscent 
of the length scale that emerged for example in 
the autocorrelation function of the non-affine field 
in recent studies of similar Lennard- Jones glasses 
|38| . These findings are in line with observations of 
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similar weak zones that grow and trigger the flip 
of neighbouring zones as depicted in fig. 1 of [H] 
confirming the validity of the picture of the dy- 
namics of slowly driven glassy materials as domi- 
nated by the accumulation and cascading of plas- 
tic events. Beyond this general mechanism a care- 
ful study of the spatio-temporal signal associated 
to the slow drive (7 = 10~^) shows different time 
scales r associated with an entire zoology of typi- 
cal sequences of plastic events. First one observes 
short lived local quadrupoles, typically visible only 
during one snapshot Tg < ILJU, and that do not 
trigger a cascade. Some local rearrangements seem 
to be locked and to survive for longer time intervals 
of the order of Tg ~ 10 — lOOLJU. In general this 
type of rearrangement triggers in its vicinity (ver- 
tically or horizontally) subsequent similar events. 
Sometimes as is the case on the forth snapshot of 
figure 2] this cascade feeds to the formation of a 
system spanning shear band. Finally one can asso- 
ciate also a timescalc to the global relaxation pro- 
cess of figure [3] which is here for the slow shear 
rate Tc = Jc/'j ~ 1000 LJU (7c corresponds here 
to the duration of the relaxation event, i.e. the to- 
tal strain associated to each downward slopes in the 
stress-strain response) and for the fastest shear rate 
Tc ~ 100 LJU. Note that these values of the typical 
duration of an entire relaxation process are in line 
with the values that one can compute from figure 
[TJ In this figure one also sees that for low enough 
shear rates 7 < 10^^ there is an intrinsic lifetime as- 
sociated to a plastic rearrangement process which 
is proven by the fact that the slope of the stress 
curves is proportional to the shear rate. For shear 
rates larger than 7 ~ 10^^ the relaxation strain be- 
comes larger than the typical strain between relax- 
ation events and therefore one can see this value of 
the shear rate as a mark of a transition to a different 
type of rheology also characterized by an important 
increase of the average yield stress as can be seen 
from figure [2l It is very striking that the avalanche 
like behavior seems to be somehow screened when 
the shear rate is increased. This result has been 
reported also elsewhere in atomic scale simulation 
[44] but also in mesoscopic yield stress models [45] 
and [15] . Only rarely studied in driven glassy mate- 
rial the growth of a cooperativity length scale near 
the glass transition is well known in the supercooled 
liquids literature (sec for example }47| ) or in simple 
lattice gas models [48] and have been interpreted in 
the framework of facilitated dynamics (for a review 
see [49]). Of course the shear that one applies on 
glassy systems breaks the symmetry of supercooled 
liquids (this is apparent for example in the existence 
of preferred orientations for the local quadrupoles 



and for the system spanning shear bands along the 
neutral axis of the external applied strain. In super- 
cooled liquids the directions of the rearrangements 
arc isotropic.) but nevertheless it appears tempt- 
ing to find, in line with supercooled liquids, a map- 
ping between the dynamics of the sheared glass and 
a simpler facilitated model. The detailed descrip- 
tion of the elementary rearranging processes that 
we propose here should help to devise reasonable 
ingredients for these models. A first approach in 
this direction was proposed by Picard et al [45] and 
we will briefly compare our results to this model at 
the end of this paper. To conclude with the descrip- 
tion of the panels H] and [5] let us mention that these 
plastic rearrangements independently of the shear 
rate emit a transverse sound wave propagating at a 
typical transverse sound speed characteristic of the 
Lennard- Jones glass (cg ~ y/ n/p ~ 3 — 4 LJU) and 
appear in general as dark regions on the snapshots 
of panel [4] and [5] This allows to introduce a new 
timescale = L/cg ^ 10 LJU, i.e. comparable to 
the life time Tg of the elementary plastic rearrange- 
ments for a system of size L = 50. We now turn 
to the detailed study of the dynamics in the driven 
systems at the atomic level. 



3 Motion of particles 

3.1 Two time correlation functions. 

The study of two time correlation function, besides 
its direct importance to quantify the characteristic 
relaxation times, enables to make comparisons with 
simple models of the rheology. In what follows, to 
examine the dynamics of the local density and the 
relaxation associated times, we compute therefore, 
on an equal foot, the self intermediate scattering 
function Fg (k, t) , 

F3(k,t) = l^cos[k.(Zir,(i))] , (1) 

i 

and the self correlation function (as in [33]) Qs{a, t) 
defined by, 

Q^(«'^)-^E^-p(-^^)' (2) 

i 

where AYi{t) = Vi{t' + t) — Yi[t') is the displace- 
ment vector. In what follows it is the spatial and 
time average of these two function that we com- 
pute and in general in what follows we replace the 
displacement vector AYi{t) by its non-affinc contri- 
bution Z\r"°(t) (or even simply by the transverse 
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displacement) defined rigorously by, 

Arrit) - r,(t' +t)~jj'^' dt"y,it")e, - r,(t) , 

(3) 

for a shear in the x direction. In practice we find 
that within a good degree of accuracy Z\r"°(<) = 
ri{t' +t) — 'jyi{t')ex — ri{t) and we use this expres- 
sion. Figure [S] represent the function Qs both in 
function of time and strain. After a slow decrease 
at small times/strains the function Qs exhibits a 
power law decay as a function of strain e, Qs{o., e) cx; 

with /? > 0.5 signaling shear induced structural 
relaxation. This power law decay is in contrast with 
the more 'classic' self intermediate scattering func- 
tion that exhibits (not shown here) a compressed 
exponential decay Fs{ky,e) oc exp(— (e/ec)^), with 
/? > 1.0. Note also that the Qs{a,€) oc decay 
is compatible with a Gaussian distribution func- 
tion of the transverse displacements P{Ay, e) and 
with a diffusive transverse mean square displace- 
ment {Ay"^). As shown in figurelHlthe convergence to 
the quasistatic curve is verified when Qs is plotted 
against strain for values of the shear rate 7 < 10~^. 

In figure[7]we have reported the relaxation times 
ti/f. - the points of intersection of the dotted line 
with the colored curves in figure[6]verifying Qsifl, ti/e 
1/e - for different shear rates and different system 
sizes. Of course the relaxation strains and relax- 
ation times are related through the simple relation 
£i/e = ^i/e * 7 and we therefore only report the re- 
laxation times for simplicity. In |501I24| similar anal- 
ysis were reported in experiments respectively on 
foams and colloids. Looking at the relaxation time 
dependence on shear rate in these studies the au- 
thors found scalings of the form ti/g cx 7"'^ with 
f ~ 0.66 in [SD| and f ~ 0.8 in [H]. Similarly 
in extensive numerical studies [5i p 52j of sheared 
Lennard- Jones glassy materials the authors have 
computed these relaxation times without explicitly 
writing to our knowledge the functional form of the 
dependence of the relaxation time on shear rate. 

Here as shown in figure [7] we find two regimes 
: for high shear rates, 7 > 10""', the structural re- 
laxation time of the sheared glass scales with the 
global shear rate as ti/^, oc 7"'^ with i' ~ 0.63, 
while for lower shear rates, 7 < 10""', the relax- 
ation functions Qs{a,e) reach a quasistatic limit 
strain limit ei/e ~ 0.04 and therefore the asso- 
ciated relaxation times scale as <i/e oc 7"'' with 
/.t ~ 1. The crossover between the quasistatic and 
shear rate dominated regimes is size dependent with 
the transition shear rate jc being lowered as the 
size of the system is increased (not shown here). 
Our data confirm the theoretically predicted 'time- 



shear superposition principle' [531154] : when time is 
scaled by ti/,, the relaxation follows a master curve 
fs{a,t/ti/e) as shown in figure [71 It is tempting to 
try and relate as in [51 > 50 > 24] the scaling exponent 
v of the structural relaxation time to the scaling 
exponent /3 that appears in the Herschel-Bulkley 
type macroscopic rheological fiow curve of the ma- 
terial where a — ay oc 7*^ (see figure [2]) . Taking, 
as is often assumed and verified numerically [52) . 
the structural relaxation time ^i/e as proportional 
to viscosity provides an expression of an effective 
stress cTeff = /xti/e7, where /i is the macroscopic 
shear modulus. Surprisingly reporting the scaling 
of the relaxation time ti/e in this expression we see 
that for high enough shear rates 7 > 10"'' the ef- 
fective stress scales with shear rate as (Te// oc 
This is in good agreement with the global mechani- 
cal response of the material and one finds a posteri- 
ori that the two coefficients /3 and v are compatible 
with the hypothesis made above and one has in- 
deed to a good approximation (3 ~ 1 ~ u for high 
shear rates. This relation breaks down for lower 
shear rates in a regime where shear banding be- 
comes the dominant relaxation mechanism, as can 
be seen for example in the panels |4] and [5] corre- 
sponding respectively to the shear rates 7 = 10~^ 
Hid 7 = 10~^. 

3.2 Mean square displacement (MSD). 

Usually to quantify the dynamics at a particle level 
one also calculates the MSD. In our two dimensional 
simulations the diffusion along the x and y direc- 
tions are not equivalent. Indeed while the diffusion 
in the shear direction (in our case the x axis) is en- 
hanced by the shear, the diffusion in the transverse 
shear-gradient direction ( y axis) is unaffected. Here 
we therefore present the MSD {Ayit)"^) in the trans- 
verse direction. Of course in the presence of rigid 
boundary conditions the diffusion along the y axis 
is limited by the presence of walls and one must be 
cautious to compute the average MSD sufficiently 
far away from the boundaries. Figure [8] shows the 
typical transverse MSD for a system containing 625 
particles under RWBCs and averaged over a total 
cumulative strain of 200% for each of 24 glass sam- 
ples. Moreover in order to avoid boundary effects 
the average is computed over one third of the sam- 
ple in the central region. Larger samples as well as 
LEBCs yield similar results and we have not re- 
produced these here for clarity. From figures [8] we 
see that MSD exhibits a transition from ballistic 
motion at short times ((Z\y^) oc t^) to diffusive mo- 
tion {{Ay"^) oc t^) for larger times. For high shear 
rate values (7 > 10^"' one can rescale the entire 
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MSD curves on a master curve git/tMSo), while 
for smaller shear rates the scaling doesn't hold for 
small times. The times tMSD are defined here as 
the intersection of the MSD curves with a 'Linde- 
mann' like criterion defined at {Aip) = 0.14 as in 
[50] . The time Imsd is a nonlinear function of shear 
rate and follows the same trend as ti/^ defined ear- 
lier (see figure |6]), but with a slightly different ex- 
ponent tMSD oc 7""^, with 1^2 ~ 0.5 (see doted line 
in the inset of figure |5]). 



3.3 DifTusion coefficient. 

The scaling of the MSD curves at different shear 
rates allows to identify the dependence of the trans- 
verse diffusion coefficient, defined by {Ay'^) = 2Dyt, 
with shear rate 7, Dy cx I/Imsd oc . In order to 
allow comparison with diffusion in quasistatically 
sheared glasses we follow Lemaitre |55] and com- 
pute the quantity Deff{Aj) = {Ay'^) /2Aj which 
is related to the usual transverse diffusion coeffi- 
cient through Deff = Dyj. In figure [S] we plot the 
effective diffusion coefficient Deff{A'^) for the var- 
ious system sizes and for various finite shear rates 
as well as under a quasistatic protocol. For all sys- 
tem sizes and shear rates we see that Deff{A'y) 
increases from a finite initial value (that increases 
with decreasing shear rate) to an asymptotic value 
for large strain. The transient strain interval ap- 
pears not strongly shear rate or size dependant and 
is of the order of ttransient = 0.25. The depen- 
dence of the asymptotic values of D^ff (we will 
call this asymptote D from now on) on shear rate 
is shown in figure IHl and displays for high shear 
rates (7 > 10"''), as expected from the relation 
Df.ff = Dy-^, the scaling D oc j"^ (z/2 ~ 0.5). 

Lemaitre [H] finds the following two limiting 
scaling behaviors of the effective diffusion coefficient 
Deff. First in the high shear rate limit one gets 
uncorrelated localized plastic events and 

Deff (X ln(L) (4) 

whilst at low shear rates Lemaitre predicts a 
linear scaling 

Deff oc L. (5) 

This linear scaling was also obtained numeri- 
cally in [Sn]. In this article the linear scaling of 
Deff is recovered by the authors with a simple 
argument if one assumes that the mechanical re- 
sponse of the material is dominated by uncorrelated 
system spanning slip lines. One then simply has 
Deff = (Ar'^/Aj^) w {A-f / {a/ L))a^/ 12, where a 



is the slip amplitude, Aj is the applied strain in- 
crement and a? 112 is the average mean square dis- 
placement associated to an individual slip line and 
accumulated during a strain of ^ a/ L. Note that a 
is assumed to be size independent. 

The difference between Maloney's and Lemaitre's 
approach is in the fact that while the later assumes 
the elementary constituents of the response to be 
the avalanches observed in deformed glassy mate- 
rials the former believes that one must take into 
account the entire system spanning slip lines that 
are formed by a cluster of avalanches as elementary 
constituents of the rheology. 

Due to the large fluctuations of the effective 
diffusion our results (figure E]) do not allow us to 
resolve clearly these questions. Therefore while we 
can not rule out the three main observations made 
by Lemaitre we can not either confirm them. First 
pertaining to the system size dependence of D at 
low shear rates we indeed observe that D grows 
with system size but our number box sizes and the 
uncertainty for each measure of D stops us from de- 
scriminating between a D^f / oc ln(L) or a D^f f L 
scaling. Second, at higher shear rates, with a sim- 
ple argument based on the long range elastic prop- 
agator of the local quadrupoles Lemaitre predicts 
Deff oc ln(L). This seems in contrast with our find- 
ings where Dgff at high shear rates seems system 
size independent (see figure |9] and compare with 
figure 5 in [44] ) . Finally due to the very long simu- 
lation time required to produce runs for shear rates 
below 7 ~ 10^^ we can not extract clearly the 
system dependence of the critical shear rate 7c sep- 
arating the system size dominated regime from the 
shear rate dominated regime. In line with previ- 
ous numerical studies (for example [44 l [57 ] [5T] ) this 
change of behavior for the three system sizes pre- 
sented here is located in the range 10^^ ^ 7 ^ 
10~^. While these results convincingly illustrated 
the influence of shear rate on the atomic motion in 
a sheared model glass they also call for extended 
simulation runs. In the next section we focus on 
what is thought to be associated with this change 
of behavior namely the existence of a growing dy- 
namical heterogeneity length scale as the shear is 
reduced. 

4 Dynamical heterogeneity 
4.1 General study 

As discussed in the introduction the dynamical het- 
erogeneity is quantified in supercooled liquids near 

^ A rmi at 7 = 10~^ for a system of size L = 100 
takes of the order of a few days for 400% strain. 



M. Tsamados: Plasticity and dynamical heterogeneity in driven glassy materials 



7 



the glass transition and more recently in jammed 
systems near the jamming point by the so-called 
four-point correlation function. Here we propose to 
extend these approaches to the case of sheared glasses 
where instead of T, the temperature, in the case 
of supercooled liquids or 0, the volume fraction, 
in granular materials we consider here 7 as the 
new control parameter. The analytical framework 
allowing to quantify the dynamical heterogeneity 
remains nevertheless identical and we therefore use 
these tools in our present analysis, in particular 
our analysis parallels the experimental study by 
Lechenault et al in j43j on the critical scalings and 
heterogeneous dynamics near the jamming/rigidity 
transition of a granular material. 

The dynamical cooperativity is quantified as the 
fluctuations of a two-point correlation function. Here 
we choose as a two point correlation function the 
transverse part self correlation function (as in equa- 
tion [2]) that we express here rather than in function 
of time in terms of the strain interval e as, 

Qsia,e)^—}_^cxpi , (6) 

i 

where a is a characteristic length scale over which 
the dynamics is probed. In figure [10] we represent 
the dependence of the spatial (index i for each parti- 
cle) and temporal/strain (index e) average Q{a, e) ~ 
((5s(a, e))i_e as a function of both the parameter a 
and the strain interval e, or rather 7 = e/2. The 
function is here plotted for a shear rate of 7 = 10^^ 
and for a sample containing 2500 particles under 
RWBCs. Q{a, e) takes values in the range [0,1], with 
Q 1 typically when the transverse motion is small 
relatively to a. Ay <ti a, and Q in the opposite 
situation when Ay ^ a. Following [43], we super- 
impose in figure [10] on the colormap of Q the root 
of the transverse MSD, ^(7,7) = y^Ay^^j). Inter- 
estingly, as for granular materials around the jam- 
ming volume fraction we sec that the MSD follows 
very well the decay of Q{a, e) and that here also the 
function Q{a, e) can be rescaled for all shear rates 
as Q{a,e) = Q'{S{€,j)/a) showing that the MSD 
defines the only microscopical relevant distance for 
a given strain and shear rate. 

Turning now to the fluctuations of the self cor- 
relation function Qs{a,e) wc define the four-point 
correlation function X4(a, e) as, 

X4(a, e)^N [{Qs{a, ef),^, - {Q,{a, e))^ . (7) 

Again we remind the reader that Xii^.!^) gives 
an estimate of the number of particles that move co- 
operatively when the sample is subjected to global 



strain e. We produce an example of this function 
for a shear rate of 7 = 10~^ in figure [TU] where 
we see that at this shear rate the maximum co- 
operativity is of the order of 5 particles. Again in 
line with [43) we obtain the same scaling with the 
mean square displacement of the four-point cor- 
relation function X4{a,e) that can be rescaled as 
X4(a, e) = /i(7, e)x4((5(e, 7)/a) where the amplitude 
/i(7, e) depends both on shear rate and strain inter- 
val. From this figure we can not determine if there 
is a point (log(a),log(e)) with finite values in this 
map corresponding to an absolute maximum of the 
function X4(o, e) as is observed by Lechenault et al. 
|43j or if the maximum is pushed at non-finite val- 
ues. 

We now turn to the influence of the shear rate on 
the dynamical cooperativity in the driven glasses. 
Figure [TT] shows for different shear rates the build 
up of cooperativity in a glass sample containing 
2500 particles. The curves start from a low value 
of X4(flj£) 9't low strain go through a maximum 
^max^^-^ at a corresponding strain e^^^{'j) {V^^^ = 
e™°^/7) and decay to zero for larger strains. Here 
the log-log representations allows to see that the 
growth with strain is of the form Xi{^) for high 
shear rate values while towards the quasistatic limit 
the behavior changes toward a Xii.^) growth. 
Note that the behavior is consistent with a bal- 
listic regime while is consistent with a regime 
dominated by collectively rearranging regions (see 
[32 )• The curves of figure [TT] allow us to extract 
the dependence of the time scale with shear 
rate 7 and we display this in the inset of figure 
[HI We find a scaling of the form oc 7""^^ 

with z/3 ~ 0.65, hence very close to the coefficient 
1^1 ^ 0.63 observed for the relaxation time ti/e as- 
sociated to the correlation functions Qs- This time 
dependence differs on the other hand slightly from 
the time Imsd extracted from the Lindemann cri- 
terion on the transverse mean square displacement 
^ 0.5 (see figure [8]). 

Again it is satisfactory to observe, in figure [TTJ 
the convergence of the X4(ajf) associated to the 
finite shear rate deformation runs towards the qua- 
sistatic data as the two sets of simulations are pro- 
duced from completely independent codes and pro- 
cedures. For the system size analyzed in this fig- 
ure (L=50) we see that the number of particles 
moving in a cooperative manner is of the order of 
^max ^ j^g jj-^ ^j^g quasistatic regime. In figure [T^ we 
have collected all the values of x™"^ for the various 
system sizes, boundary conditions and shear rates 
that we have analyzed. The result shows strikingly 
the growth of cooperative length scale with decreas- 
ing shear rate for all system sizes. This plot illus- 
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trates again (as was the case for the flow curve or 
for the diffusion properties) two regimes, namely a 
high shear rate regime (above a critical shear rate 
7 ^ 7c) where xT'^^ grows with decreasing 7 and 
a plateau regime where xT"'^ saturates to a system 
size limited value. While the data are still quite 
scattered they allow to extract a typical scaling 
coefficient ^ for the dependence of xT°'^ on. 7 as 
^moa: ^ ^-^^ with ^ ^ OA — 0.6. The dependence 
of x™"^ on system size in the quasistatic regime is 
too noisy to be quantified precisely at the present 
stage of our study. It is obvious nevertheless that 
the dynamical cooperativity x™°^ grows in the qua- 
sistatic regime with system size in a way that in- 
dicates finite size effects. We will come back to the 
scaling exponent fi in the last section of this paper 
and relate our findings with previous observations. 
But let us try first to visualize the spatial structures 
associated to this build up of a growing dynamical 
heterogeneity length scale as 7 — )■ 0. 

4.2 Spatio-temporal structures and a simple 
model for the cooperativity length 

There are two aspects that we would like to il- 
lustrate here. First how does the dynamical het- 
erogeneity build up during the strain of a glassy 
material, in other words when one looks at spatial 

maps of Ql = exp (^ '^^2!^^ ) ^or increasing values 
of e. Second how the dynamical heterogeneity is af- 
fected by the value of the imposed shear rate, in 
other words when one plots spatial maps of Ql at 
^civu t)ut for different shear rates. The next two 
panels of fT3l and [Ml illustrate respectively the spa- 
tial fluctuations of exp (^ ^2aif^ ) ^^'^ ^-^P ( ^2a^^ ) 
for various strain intervals, while the panel [15] illus- 
trates the fluctuations of exp f^^^^ at the peak 
of the four-point correlation function Xi for differ- 
ent shear rates. Comparing figure [T3l and figure [T4l 
confirms the already observed fact that the relax- 
ation of the Ql to zero is faster in the x direction 
that in the y direction. This anisotropy is here en- 
hanced by the presence of walls but would be visi- 
ble also under LEBCs. It is striking to sec in figure 
[T5]the growth of a cooperativity length scale as the 
shear rate is decrease from 10""^ to I0~^. Indeed one 
sees that the response of the glass to the external 
shear rates becomes more and more heterogeneous 
as the shear rate is lowered and that at 7 = 10"^, 
for example, the typical size of the clusters of parti- 
cles that have moved more than a = 0.1 (the white 
particles) seems to reach an important fraction of 
system size. Note that in all these maps the clusters 



seem to form anisotropic structures elongated along 
the y axis. We attribute this apparent anisotropy to 
the formation of vortices. 

Let us finally conclude this article by giving a 
simple physical argument that allows to understand 
the scaling of this cooperative length scale with 
shear rate. The argument goes as follows. Assume 
that the regions that are prone to fail plastically (we 
call them the triggering points) under shear are ho- 
mogeneously distributed in the glass and that the 
density per unit strain pte (te stands for triggering 
events) of these points is a constant that is inde- 
pendent of strain rate. Then during a time interval 
t there are 7 x t x pte triggering points that arc 
excited and one can define the average distance Ite 
between these points to be Ite = (pteji) ^^'^ , where 
d is the dimension. These points by definition are 
triggering a quadrupolar event that as we have dis- 
cussed in the previous section can induce further 
plastic rearrangements in its vicinity. More specifi- 
cally a plastic rearrangement induces a quadrupo- 
lar redistribution of the stress in its surrounding 
and one expects an increased probability of having 
a new plastic event where the stress is increased, 
i.e. along the vertical and horizontal axis. The sim- 
plest hypothesis concerning the propagation of the 
events is that it occurs through a diffusive proccsfl 
One can then define a new length scale Id = V Dt, 
where D is the diffusion coefficient. Now to extract 
a characteristic length scale we identify these two 
length scales, he = l-D which yields a characteristic 
time tc, 

-1 -1 -1 

^ pd/2 + 1 JJl + 2/djd/2+l ^ (8) 

and a corresponding length scale scaling as, 

Ic = p"^ D^'j'^ . (9) 

In two dimensions this corresponds to the scalings 
tc (X 7~^/^ and Ic cx; 7"^/^ while in three dimensions 
one expects te oc 7^^/^ and le oc 7^-'^/^. These re- 
sults argue rather well with the observations made 
earlier on the relaxation time scale and on the growth 
of a cooperativity length scale. We argue that this 
length should correspond to the maximum extent of 
the plastic cascade and therefore to the maximum 

^ One could make the model more quantitative by 
mapping it to the problem of diffusion on a grid. One 
can estimate the diffusion coefficient from atomic con- 
siderations as D C^/''') where ^ is the optimal distance 
between successive plastic rearrangements and r is the 
duration of plastic event. From figure|4]one can estimate 
^ to be of the order of ^ ~ 20 — 30cr and one can identify 
T with the duration of a plastic event r ~ ~ 10 — 100. 
This gives 1 < -D < 100 
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cooperativity length. Indeed when Ijj becomes of 
the order of the distance between triggering events 
ke each plastic cascade starts to 'feel' the presence 
of the neighboring cascades and its progression is 
perturbed. Finally wc sec that the exponent —1/4 
is in good agreement with what is observed numer- 
ically in figure [T^] but also with the predictions of 
the kinetic clastoplastic model introduced in [581 
as well as with the prediction of mesoscopic elasto- 
platic models jl5] . 

5 Conclusion 

The rheology of a model two-dimensional glass was 
discussed with an emphasis on the complex plastic 
response and dynamical heterogeneity that builds 
up in the material under shear. The detailed analy- 
sis of the plastic deformation at different shear-rates 
showed that the glass follows different flow regimes. 
Below a system size dependent critical shear-rate 
the mechanical response reaches a quasistatic limit 
characterized by finite size effects, cascades of plas- 
tic rearrangements and yield stress, while at higher 
shear rates the rheological properties are determined 
by the externally applied shear-rate. In the later 
regime wc reported on the growth of a coopera- 
tivity length scale and discussed the scaling of this 
length with shear-rate. We proposed a simple model 
based on the diffusive propagation of plastic events 
that reproduces this scaling. The glassy system was 
shown to develop a diverging cooperativity length 
scale with lowering shear rate and one can ask how 
this dynamical length scale affects the fiow of con- 
fined glassy systems. Preliminary results concern- 
ing the effect of temperature have shown surpris- 
ing effects associated with the presence of a small 
but finite temperature and many further simula- 
tion runs at various temperatures and shear rates 
arc needed to apprehend quantitatively the rheol- 
ogy of the glasses. Similarly our results on a two 
dimensional system call for a generalization to the 
three dimensional case. Our results have also con- 
firmed at the atomic level many of the predictions of 
the mesoscopic extremal elasto-plastic models [59| 
145] and quantitative comparison would require an 
extended amount of simulated material, larger sys- 
tem sizes and lower shear rates. Finally our work 
is a first step towards a better understanding of 
the elementary building blocks needed to construct 
a mesoscopic description of the rheology of glassy 
materials along with innovative constitutive laws. 
Another possible extension of this work would be 
to map the rheology of the structural glasses ob- 
tained through the use of atomic scale simulations 
to a kinetically constrained model, the collaborative 



nature of the dynamical building blocks presented 
in this article seems indeed to validate this type of 
approach. 

During the course of this work, I have benefited 
from discussions with Glaus Heussinger, Pinaki Chaud- 
huri, Lyderic Bocquet, Anne Tanguy and Jean-Louis 
Barrat. Numerical computations were carried out 
on P2CIIPD (University Lyon) computers. 
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Fig. 1. Left : stress strain mechanical response of a glass sample containing 625 particles sheared under RWBCs 
at shear rates ranging from 7 — 10~^ to 7 = 10~^ (thin colored lines from top to bottom, 7 — 10^^ (black), 
7 = 5- 10"^ (red), 7 = 2.5 • 10"^ (green), 7 = 10~^ (dark green), 7 = 5- 10"* (brown), 7 = 2.5 • 10"'' (gray), 
7 = lO""* (violet), 7 = 5- lO"'^ (cyan), 7 = 2.5 ■ 10"'' (magenta) and 7 = lO"'' (orange)). The thick black 
line corresponds to the quasistatic shear protocol. Right : Zoom in a portion of the total mechanical response, 
illustrating the typical relaxation time associated to a plastic rearrangement in the glass signaled in the quasistatic 
protocol by an abrupt stress drop. 
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shear rate Sheai^ rate 

Fig. 2. Flow curves associated with the different system sizes. Small open symbols correspond to finite shear 
rate values ranging from 10~^ to 10~^. The larger solid points on the vertical axis correspond to the quasistatic 
protocol 7 = 0. Left : Shear under LEBCs with fit parameter to the Herscel-Bulkley rheological law, ty — 0.32, 
ci = 7.2 and /3 = 0.38. Right : Shear under RWBCs, with the dashed line representing a Herschel-Bulkley fit 
T = TY + ci-yl^, with TY ~ 0.36, ci = 7.4 and /? = 0.4. 
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Fig. 6. Left : 2-point correlation function Qs{a,e) — '}2i^'^Pi~ 2a'^ ) against strain e for a 
system with LEBC. Thin lines (color on pdf version) : from top to bottom shear rate values are 
0.01,0.005,0.0025,0.001,0.0005,0.00025,0.0001,0.00005 and 0.000025. Dashed line : Qs{a,e) vs e for quasistatic 
shear. Doted line corresponds to the value 1/e for which we calculate the 1/e-relaxation strains and times. Right: 
Qs{a,t) against time t. In both figures the parameter a is chosen equal to 0.1. 




Fig. 7. Left : Relaxation time fi/^ vs shear rate 7 for the system of figure[6l Dotted line t-i/^ oc 7~^, dashed line 
oc -y""'^"^. Right: Rescaled self correlation function Qs{a,t) of figure|B]when time is rescaled by the structural 
relaxation time ti^^. All curves superimpose rather well on a master curve fs{a,t/ti^^). 
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Fig. 8. Left : Transverse mean square displacement (MSD) {Ay^) versus time for different shear rates 7 (same 
color code as in previous figures). The dashed line marks the arbitrarily chosen criterion distance a verifying 

= 0.14 corresponding to a 'Lindemann' criterion for the transverse MSD of the particles. The intersection of 
this line with the colored curves marks the times Imsd- The two thick black lines correspond to the power laws 

and , i.e. respectively to the ballistic and diffusive regime. Right: same figure where time is rescaled by the 
times Imsd- These plots are obtained in configurations with RWBCs and the average are computed in a central 
region of the samples therefore avoiding direct influence of the walls. 
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Fig. 9. Left : Efltective diffusion coefficients defined as D^ffiAj) = (Ay^) /2A'y for different shear rates ranging 
from bottom to top from -y = 10~* to 10~^ (same color code as previous figures) and for the quasistatic shear 
protocol (black dashed line). De// is computed here for sample containing 2500 particles under RWBCs. The 
spatial averaging is performed sufficiently far from the rigid walls to avoid boundary effects. Right: The asymptotic 
values Def / are plotted for various shear rates, system sizes and boundary conditions. The dashed line marks the 
power law De// oc j'"'^ as a guide to the eyes. 




Fig. 10. Dynamical correlation functions computed over the particles of sample containing 2500 particles and 
sheared at 7 = 10"'' under RWBCs for a total strain etot ~ 200%. Left : Correlation function Qs (a, 7) as a function 
of the probing length a and the strain 7 in a log-log colormap. Right: Four-point correlation function Xii^-jt) 
a log-log colormap. 
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Fig. 11. Four-point correlation function as a function of strain e for different shear rates (same color code as 
previous figures) computed on samples containing 2500 particles and sheared under RWBCs. The averaging is 
done over 15 samples over a total strain of e = 200% on each configuration. Insert : Shear rate dependence of the 
time t™^ associated to the peak of X4- The different symbols correspond to different system sizes and boundary 
conditions, see figure |8] for the legend. The dashed line corresponds to the power law 7"°'^ and is shown as a 
guide to the eyes. 
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Fig. 12. maximum values x?""^ of the four-point correlation function as a function of shear rate 7 for various 
system sizes and boundary conditions. The dashed line shows the power law oc 7~'^, with fi = 0.6. 
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M. Tsamados: Plasticity and dynamical heterogeneity in driven glassy materials 





Fig. 13. Top left corner : X'ti.^)- The red marks correspond to the strain intervals at which the five spatial 

maps of the self correlation function Ql{e) = exp ) computed. From top to bottom and from left to 

right, £ = 10~^,£ = 2 ■ 10~^,e — 4 ■ 10~^,e = 8 • 10~'^,e = 16 ■ 10~^. All figures are obtained on a sample containing 
10000 particles under RWBCs and at a shear rate of 7 = 5 ■ lO"**. 



M. Tsamados: Plasticity and dynamical heterogeneity in driven glassy materials 



21 






Fig. 14. Top left corner : Xi{^)- The red marks correspond to the strain intervals at which the five spatial maps 

of the self correlation function Ql{e) — exp ^ ^ 2'al'' ^ ^-re computed. From top to bottom and from left to right, 

e = 2.5 • 10~^,e = 2.5 • 10~^,e = 5 • 10~^,e = 10"^, e = 2 • 10^^. All figures are obtained on a sample containing 
10000 particles under RWBCs and at a shear rate of 7 = 5 ■ 10^''. 
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Fig. 15. Top left corner : X'^i^)- The red marks correspond to the strain intervals at which the five spatial 

maps of the self correlation function Ql{e) = exp ^ ^2al^ ) computed. From top to bottom and from left to 

right, 7 = 10~^,7 = 5 ■ 10^'*, 7 = 10^'*,7 — 5 • 10~^,7 = 10~^. All figures are obtained on a sample containing 
10000 particles and for RWBCs. 
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